/*
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */

package jmeshless;


import java.io.FileInputStream;
import java.io.FileNotFoundException;
import java.util.ArrayList;
import java.util.Scanner;

import org.jfree.data.xy.XYSeries;
import org.jfree.data.xy.XYSeriesCollection;
import org.yourorghere.Node;
import org.yourorghere.Triangle;

import Integration.IFunction2D;
import MF1DSolver.MF1DSolver;
import MeshlessShapeFunctions.JMeshlessInterpolationFunction1D;
import MeshlessShapeFunctions.JMeshlessInterpolationFunctionLiu;
import MeshlessShapeFunctions.JStKernel;
import MeshlessShapeFunctions.Particle;

/**
 *
 * @author mark
 */
public class Main {

    /**
     * @param args the command line arguments
     */
    public static void main(String[] args) throws FileNotFoundException {

        /*ArrayList<Particle> ls = new ArrayList<Particle>();
        ls.add(new Particle(0.0, 0.0, 1));
        ls.add(new Particle(0.0, 0.25, 1));
        ls.add(new Particle(0.0, 0.75, 1));
        ls.add(new Particle(0.0, 1.0, 1));

        ls.add(new Particle(0.0, -0.25, 1));
        ls.add(new Particle(0.0, -0.75, 1));
        ls.add(new Particle(0.0, -1.0, 1));

        ls.add(new Particle(0.25, 0.0, 1));
        ls.add(new Particle(0.25, 0.25, 1));
        ls.add(new Particle(0.25, 0.75, 1));
        ls.add(new Particle(0.25, 1.0, 1));

        ls.add(new Particle(0.25, -0.25, 1));
        ls.add(new Particle(0.25, -0.75, 1));
        ls.add(new Particle(0.25, -1.0, 1));

        ls.add(new Particle(-0.25, 0.0, 1));
        ls.add(new Particle(-0.25, 0.25, 1));
        ls.add(new Particle(-0.25, 0.75, 1));
        ls.add(new Particle(-0.25, 1.0, 1));

        ls.add(new Particle(-0.25, -0.25, 1));
        ls.add(new Particle(-0.25, -0.75, 1));
        ls.add(new Particle(-0.25, -1.0, 1));


        ls.add(new Particle(0.75, 0.0, 1));
        ls.add(new Particle(0.75, 0.25, 1));
        ls.add(new Particle(0.75, 0.75, 1));
        ls.add(new Particle(0.75, 1.0, 1));

        ls.add(new Particle(0.75, -0.25, 1));
        ls.add(new Particle(0.75, -0.75, 1));
        ls.add(new Particle(0.75, -1.0, 1));

        ls.add(new Particle(-0.75, 0.0, 1));
        ls.add(new Particle(-0.75, 0.25, 1));
        ls.add(new Particle(-0.75, 0.75, 1));
        ls.add(new Particle(-0.75, 1.0, 1));

        ls.add(new Particle(-0.75, -0.25, 1));
        ls.add(new Particle(-0.75, -0.75, 1));
        ls.add(new Particle(-0.75, -1.0, 1));

        ls.add(new Particle(1.0, 0.0, 1));
        ls.add(new Particle(1.0, 0.25, 1));
        ls.add(new Particle(1.0, 0.75, 1));
        ls.add(new Particle(1.0, 1.0, 1));

        ls.add(new Particle(1.0, -0.25, 1));
        ls.add(new Particle(1.0, -0.75, 1));
        ls.add(new Particle(1.0, -1.0, 1));

        ls.add(new Particle(-1.0, 0.0, 1));
        ls.add(new Particle(-1.0, 0.25, 1));
        ls.add(new Particle(-1.0, 0.75, 1));
        ls.add(new Particle(-1.0, 1.0, 1));

        ls.add(new Particle(-1.0, -0.25, 1));
        ls.add(new Particle(-1.0, -0.75, 1));
        ls.add(new Particle(-1.0, -1.0, 1));*/
        ArrayList<Particle> ls1D = new ArrayList<Particle>();
        ArrayList<JMeshlessInterpolationFunction1D> MeshlessShapeList = new ArrayList<JMeshlessInterpolationFunction1D>();
        
//        ls1D.add(new Particle(0.0, 0.0, 0));
//        h = 0.133;
//        MeshlessShapeList.add(new JMeshlessInterpolationFunction1D(1, ls1D.get(ls1D.size()-1).getX(), 2.33*h, new JGaussian()));
//        ls1D.add(new Particle(0.1, 0.0, 1));
//        MeshlessShapeList.add(new JMeshlessInterpolationFunction1D(1, ls1D.get(ls1D.size()-1).getX(), 2.33*h, new JGaussian()));
//        ls1D.add(new Particle(0.2, 0.0, 2));
//        MeshlessShapeList.add(new JMeshlessInterpolationFunction1D(1, ls1D.get(ls1D.size()-1).getX(), 2.33*h, new JGaussian()));
//        ls1D.add(new Particle(0.3, 0.0, 3));
//        MeshlessShapeList.add(new JMeshlessInterpolationFunction1D(1, ls1D.get(ls1D.size()-1).getX(), 2.33*h, new JGaussian()));
//        ls1D.add(new Particle(0.4, 0.0, 4));
//        MeshlessShapeList.add(new JMeshlessInterpolationFunction1D(1, ls1D.get(ls1D.size()-1).getX(), 2.33*h, new JGaussian()));
//        ls1D.add(new Particle(0.5, 0.0, 5));
//        MeshlessShapeList.add(new JMeshlessInterpolationFunction1D(1, ls1D.get(ls1D.size()-1).getX(), 2.33*h, new JGaussian()));
//        ls1D.add(new Particle(0.6, 0.0, 6));
//        MeshlessShapeList.add(new JMeshlessInterpolationFunction1D(1, ls1D.get(ls1D.size()-1).getX(), 2.33*h, new JGaussian()));
//        ls1D.add(new Particle(0.7, 0.0, 7));
//        MeshlessShapeList.add(new JMeshlessInterpolationFunction1D(1, ls1D.get(ls1D.size()-1).getX(), 2.33*h, new JGaussian()));
//        ls1D.add(new Particle(0.8, 0.0, 8));
//        MeshlessShapeList.add(new JMeshlessInterpolationFunction1D(1, ls1D.get(ls1D.size()-1).getX(), 2.33*h, new JGaussian()));
//        h = 0.07;
//        ls1D.add(new Particle(0.85, 0.0, 9));
//        MeshlessShapeList.add(new JMeshlessInterpolationFunction1D(1, ls1D.get(ls1D.size()-1).getX(), 2.33*h, new JGaussian()));
//        ls1D.add(new Particle(0.9, 0.0, 10));
//        MeshlessShapeList.add(new JMeshlessInterpolationFunction1D(1, ls1D.get(ls1D.size()-1).getX(), 2.33*h, new JGaussian()));
//        h = 0.035;
//        ls1D.add(new Particle(0.925, 0.0, 11));
//        MeshlessShapeList.add(new JMeshlessInterpolationFunction1D(1, ls1D.get(ls1D.size()-1).getX(), 2.33*h, new JGaussian()));
//        ls1D.add(new Particle(0.95, 0.0, 12));
//        MeshlessShapeList.add(new JMeshlessInterpolationFunction1D(1, ls1D.get(ls1D.size()-1).getX(), 2.33*h, new JGaussian()));
//        ls1D.add(new Particle(0.975, 0.0, 13));
//        MeshlessShapeList.add(new JMeshlessInterpolationFunction1D(1, ls1D.get(ls1D.size()-1).getX(), 2.33*h, new JGaussian()));
//        ls1D.add(new Particle(1.0, 0.0, 14));
//        MeshlessShapeList.add(new JMeshlessInterpolationFunction1D(1, ls1D.get(ls1D.size()-1).getX(), 2.33*h, new JGaussian()));

//        for (int i = 0; i < ls1D.size()-1; i++) {
//            h = -ls1D.get(i).getX() + ls1D.get(i+1).getX();
//            MeshlessShapeList.add(new JMeshlessInterpolationFunction1D(1, ls1D.get(i).getX(), 1.33*h, new JGaussian()));
//        }
//        MeshlessShapeList.add(new JMeshlessInterpolationFunction1D(1, ls1D.get(ls1D.size()-1).getX(), 1.33*h, new JGaussian()));

        int nn = 25;
        double a =0,b=1;
        double h = (b-a)/(double)nn;

        for (int i = 0; i <=nn; i++) {
            ls1D.add(new Particle(a+h*i, 0, i));
            ls1D.get(ls1D.size()-1).setRho(2.25*h);
            //MeshlessShapeList.add(new JMeshlessInterpolationFunction1D(1, a+i*h, 1.33*h, new JExponential(0.3)));
            MeshlessShapeList.add(new JMeshlessInterpolationFunction1D(1, a+i*h, 2.25*h, new JStKernel(2.25*h)));
            //MeshlessShapeList.add(new JMeshlessInterpolationFunction1D(1, a+i*h, 1.33*h, new JStKernel()));
//            System.out.println(ls1D.get(i).getX()+"\t");
        }

        System.out.println();
        System.out.println("=========================================================");

        MF1DSolver solver = new MF1DSolver();
        solver.setLeft(0.0);
        solver.setRight(1.0);
            solver.setfString("100");
            solver.setwString("10");
            solver.setqString("0");
            solver.setpString("1");
        solver.setListOfShapeFunction(MeshlessShapeList);
        solver.setNodeList(ls1D);
        solver.generate();
        solver.setAVal(0.0);
        solver.setDirihetA(true);
        solver.setBVal(0.0);
        solver.setDirihetB(true);
        double[] res = solver.solvePenalty();
//        double[] res = solver.solveLagrangeMulti();
        double x = 0;
        while(x<=1.0001){
            double tmp = 0.0;
            for (int j = 0; j <MeshlessShapeList.size(); j++) {
                tmp+= MeshlessShapeList.get(j).calculate(x)*res[j];
            }
//            System.out.print(x+"\t"+tmp+"\n");
            x+=0.01;
        }
        
        
        double tmp = 0.0;
        x = 0.0;
        for (int j = 0; j <MeshlessShapeList.size(); j++) {
            tmp+= MeshlessShapeList.get(j).calculate(x)*res[j];
        }
//        System.out.print(x+"\t"+tmp+"\n");
        System.out.printf("%4.8e", tmp);
        System.out.println();
        x = 1.0;
        tmp = 0.0;
        for (int j = 0; j <MeshlessShapeList.size(); j++) {
            tmp+= MeshlessShapeList.get(j).calculate(x)*res[j];
        }
//        System.out.print(x+"\t"+tmp+"\n");
        System.out.printf("%4.8e", tmp);
        XYSeries series1 = new XYSeries("розвязок методом штрафу", true, true);
        x = 0;
        while(x<=1.001){
            tmp = 0.0;
            for (int j = 0; j <MeshlessShapeList.size(); j++) {
                tmp+= MeshlessShapeList.get(j).calculate(x)*res[j];
            }
            series1.add(x,tmp);
            x+=0.0001;
        }

        NewOkCancelDialog dialog = new NewOkCancelDialog(null, false, new XYSeriesCollection(series1));
        dialog.setLocationRelativeTo(null);
        //dialog.setUndecorated(true);
        dialog.setVisible(true);
        
        
        
        //for (int i = 0; i <=nn; i++) {
            //System.out.print(ls1D.get(i).getX()+"\t");
        //    double tmp=0;
        //        for (int j = 0; j < MeshlessShapeList.size(); j++) {
        //        JMeshlessInterpolationFunction1D jMeshlessInterpolationFunction1D = MeshlessShapeList.get(j);
        //        tmp+=jMeshlessInterpolationFunction1D.calculate(ls1D.get(i).getX());
        //    }
        //    System.out.print(tmp+"\t");
        //}
        //System.out.println();
        //System.out.println("=========================================================");


        //for (int i = 0; i <=nn; i++) {
            //System.out.print(ls1D.get(i).getX()+"\t");
        //    JMeshlessInterpolationFunction1D jMeshlessInterpolationFunction1D = MeshlessShapeList.get(i);

        //    System.out.print(jMeshlessInterpolationFunction1D.calculate(jMeshlessInterpolationFunction1D.getXii())+"\t");
        //}

        //jm1D = new JMeshlessInterpolationFunction1D(6, 0.0, 1.0, new JGaussian());
        //JMeshlessInterpolationFunction1D2 Mes = new JMeshlessInterpolationFunction1D2(6, -1.5, 1.0, new JGaussian());
        //jm1D.setNodeList(ls1D);
        //Mes.setNodeList(ls1D);

        /*nn = 250;
        h = (b-a)/(double)nn;
        ls = new ArrayList<Particle>();
        for (int i = 0; i <= nn; i++) {
        ls.add(new Particle(a+h*i, 0, i));
        System.out.print(ls.get(i).getX()+"\t");
        }*/
        System.out.println();
        System.out.println("=========================================================");
        //IKernelFunction func = new JExponential(0.3);
        //for (int i = 0; i <= nn; i++) {
        //    System.out.print(jm1D.calculateDiffuseDerivaty(ls.get(i).getX())+"\t");
        //System.out.print(jm1D.calculate(ls.get(i).getX())+"\t");
        //System.out.print(Mes.calculate(ls.get(i).getX())+"\t");
        //System.out.print(jm1D.calculate(ls.get(i).getX())+"\t");
        //System.out.print(func.evaluate(ls.get(i).getX())+"\t");
        //}
        //System.out.println();
        //System.out.println("=========================================================");
        //System.out.println();
        //System.out.println(new GaussQuadrature().Gauss(0, 1, 0.01, new Func()));
        //System.out.println(new GaussQuadrature2D().Gauss(0, 1,0,1, 0.01, new Func2D()));
        //System.out.println(trap2D(0, 1, 0, 1, 1000, new Func2D()));
            //FEMShapeFunction1 fem = new FEMShapeFunction1(new MathLib.triangulation.SaveStructures.Node(0, 0, true, 1), new MathLib.triangulation.SaveStructures.Node(1, 0, true, 1), new MathLib.triangulation.SaveStructures.Node(0, 1, true, 1));
            //fem.setNodeNum(1);
        //System.out.println(new GaussQuadrature2D5Point().GaussOnTriangle(new MathLib.triangulation.SaveStructures.Node(0, 0, true, 1), new MathLib.triangulation.SaveStructures.Node(1, 0, true, 1), new MathLib.triangulation.SaveStructures.Node(0, 1, true, 1), 0.01, new  Func2D()));
        //System.out.println(new GaussQuadrature2D5Point().g(0,1,0,1,new  Func2D(),new MathLib.triangulation.SaveStructures.Node(0, 0, true, 1), new MathLib.triangulation.SaveStructures.Node(1, 0, true, 1), new MathLib.triangulation.SaveStructures.Node(0, 1, true, 1)));
            //System.out.println(new Gauss2D(20).gaussOnTriangle(fem,new MathLib.triangulation.SaveStructures.Node(0, 0, true, 1), new MathLib.triangulation.SaveStructures.Node(1, 0, true, 1), new MathLib.triangulation.SaveStructures.Node(0, 1, true, 1)));
        //System.out.println(new GaussQuadrature2D().GaussOnTriangle(new MathLib.triangulation.SaveStructures.Node(0, 0, true, 1), new MathLib.triangulation.SaveStructures.Node(1, 0, true, 1), new MathLib.triangulation.SaveStructures.Node(0, 1, true, 1), 0.01, fem));
        /*System.out.println();
        System.out.println("1=========================================================");

        FEMShapeFunction2 fem2 = new FEMShapeFunction2(new MathLib.triangulation.SaveStructures.Node(0, 0, true, 1), new MathLib.triangulation.SaveStructures.Node(0.5, 0, true, 1), new MathLib.triangulation.SaveStructures.Node(1, 0, true, 1), new MathLib.triangulation.SaveStructures.Node(0.5, 0.5, true, 1), new MathLib.triangulation.SaveStructures.Node(0, 1, true, 1), new  MathLib.triangulation.SaveStructures.Node(0, 0.5, true, 1));

        for(int i=1;i<=6;i++){
        for(int j=1;j<=6;j++){
        System.out.print(new Gauss2D(20).gaussOnTriangle(fem2.calculateMultFuncs(i, j),new MathLib.triangulation.SaveStructures.Node(0, 0, true, 1), new MathLib.triangulation.SaveStructures.Node(1, 0, true, 1), new MathLib.triangulation.SaveStructures.Node(0, 1, true, 1))+"\t");
        }
        System.out.println();
        }
        System.out.println("2=========================================================");
        //double[][] mas = GenerateMe(new MathLib.triangulation.SaveStructures.Node(0, 0, true, 1), new MathLib.triangulation.SaveStructures.Node(1, 0, true, 1), new MathLib.triangulation.SaveStructures.Node(0, 1, true, 1));

        Func2D ftmp = new  Func2D();
        for(int i=0;i<=5;i++){
        ftmp.i = i+1;
        for(int j=0;j<=5;j++){
        ftmp.j = j+1;
        System.out.print(new Gauss2D(20).gaussOnTriangleTMP(ftmp,new MathLib.triangulation.SaveStructures.Node(0, 0, true, 1), new MathLib.triangulation.SaveStructures.Node(1, 0, true, 1), new MathLib.triangulation.SaveStructures.Node(0, 1, true, 1))+"\t");
        }
        System.out.println();
        }
        System.out.println("3=========================================================");
        double[][] mas = GenerateMe(new MathLib.triangulation.SaveStructures.Node(0, 0, true, 1), new MathLib.triangulation.SaveStructures.Node(1, 0, true, 1), new MathLib.triangulation.SaveStructures.Node(0, 1, true, 1));
        for(int i=0;i<=5;i++){
        for(int j=0;j<=5;j++){
        System.out.print(mas[i][j]+"\t");
        }
        System.out.println();
        }*/
        //System.out.println(new GaussQuadrature(10).GaussOnGamma(new MathLib.triangulation.SaveStructures.Node(0, 0, true, 1), new MathLib.triangulation.SaveStructures.Node(0, 1, true, 1), fem) );

        //jm = new JMeshlessInterpolationFunctionLiu(1, 0.0, 0.0, 2.0, new JGaussian());
        //jm.setNodeList(ls);

        String path = "C:\\Users\\mark_o\\Documents\\NetBeansProjects\\MySimpleFEM\\New Folder\\New Folder\\";
        //getNodes(path+"\\file.1.node");
        //getTriangles(path+"\\file.1.ele");
        //getSolvedResults(path+"\\file.txt");
        //for (int i = 0; i < Results.size(); i++) {
        //    System.out.println(NodeList.get(i).x+" "+NodeList.get(i).y + " " + Results.get(i));
        //}

    }

    public static void getNodes(String FilePath) throws FileNotFoundException{
        /*Scanner in = new Scanner((new FileInputStream(FilePath)));
        NodeList  = new ArrayList<Node>();
        int n = in.nextInt();
        in.nextInt();
        in.nextInt();
        in.nextInt();
        for (int i = 0; i <n; i++) {
        in.next();
        NodeList.add(new Node( ( Double.parseDouble(in.next()) ), ( Double.parseDouble(in.next()) ) , ( isTrue(Integer.parseInt(in.next())) )  ));
        }*/
        NodeList = new ArrayList<Node>();
        double a= -3;
        double b= 3;
        double h = (b-a)/(double)100;
        for (int i = 0; i < 100; i++) {
            for (int j = 0; j < 100; j++) {
                NodeList.add(new Node(a+i*h, a+j*h, true));
            }
        }
    }

    private static boolean isTrue(int s){
        if(s==0){
            return false;
        }else{
            return true;
        }
    }

    public static void getTriangles(String FilePath) throws FileNotFoundException{
            Scanner in = new Scanner((new FileInputStream(FilePath)));
            TrList  = new ArrayList<Triangle>();
            int n = in.nextInt();
            in.nextInt();
            in.nextInt();
        for (int i = 0; i <n; i++) {
            in.next();
            TrList.add(new Triangle( ( Integer.parseInt(in.next()) ), ( Integer.parseInt(in.next()) ) , ( Integer.parseInt(in.next()) )  ));
        }
    }

    public static void getSolvedResults(String FilePath) throws FileNotFoundException{
        Results = new ArrayList<Double>();
        for (int i = 0; i < NodeList.size(); i++) {
            Node node = NodeList.get(i);
            //Results.add(jm.calculate(node.x, node.y));
            //Results.add(jm.calculateDerivatyX(node.x, node.y));
            Results.add(jm.calculateDiffuseDerivatyX(node.x, node.y));

            //FEMShapeFunction1 fem = new FEMShapeFunction1(new MathLib.triangulation.SaveStructures.Node(0, 0, true, 1), new MathLib.triangulation.SaveStructures.Node(1, 0, true, 1), new MathLib.triangulation.SaveStructures.Node(0, 1, true, 1));
            //fem.setNodeNum(2);
            //Results.add(fem.calculate(node.x, node.y));
            //System.out.println(Results.get(Results.size()-1));
        }

    }


    public static double trap2D(double a, double b,double ay, double by, int n, IFunction2D func){
            double res = 0.0;
            double h = (b - a) / (double)n;
            double hy = (by - ay) / (double)n;
            if (n > 1) {

                for (int j = 1; j < n; j++) {
                    for (int i = 1; i < n; i++) {
                        res += func.calculate(a + h * i,ay + hy * j);
                    }
                }
            }
            res += func.calculate(a,ay) / 2.0 + func.calculate(b,by) / 2.0;
            res *= h*hy;

            //System.out.println(res);
            return res;
        }


    public static double[][] GenerateMe(MathLib.triangulation.SaveStructures.Node A,MathLib.triangulation.SaveStructures.Node B,MathLib.triangulation.SaveStructures.Node C){
        double d = (2*ComputeTriangleArea(A, B, C))/72.0;
        return new double[][]{
            { 12*d,    0, -2*d, -8*d, -2*d,    0},
            { 0,    64*d,    0, 32*d, -8*d, 32*d},
            { -2*d,    0, 12*d,    0, -2*d, -8*d},
            { -8*d, 32*d,    0, 64*d,    0, 32*d},
            { -2*d, -8*d, -2*d,    0, 12*d,    0},
            { 0,    32*d, -8*d, 32*d,    0, 64*d}
        };
    }

    public static double ComputeTriangleArea(MathLib.triangulation.SaveStructures.Node A,MathLib.triangulation.SaveStructures.Node B,MathLib.triangulation.SaveStructures.Node C){

        double a =Math.sqrt( (A.x - B.x)*(A.x - B.x) + (A.y - B.y)*(A.y - B.y) );   //a-b
        double b =Math.sqrt( (C.x - B.x)*(C.x - B.x) + (C.y - B.y)*(C.y - B.y) );
        double c =Math.sqrt( (A.x - C.x)*(A.x - C.x) + (A.y - C.y)*(A.y - C.y) );
        double p = (a+b+c)/2.0;
        return Math.sqrt( p*(p-a)*(p-b)*(p-c) );
    }
    
    public static ArrayList<Double> Results;
    public static ArrayList<Node> NodeList;
    public static ArrayList<Triangle> TrList;
    public static JMeshlessInterpolationFunctionLiu jm;
    public static JMeshlessInterpolationFunction1D jm1D;
}
